#include <fstream>
#include <iostream>
#include <vector>
#include <cmath>
#include <ctime>

const int IMAGE_WIDTH = 1280;
const int IMAGE_HEIGHT = 720;

struct Pixel {
    Pixel() {}
    Pixel(float r, float g, float b, float a) : r(r), g(g), b(b), a(a) {}
    float r, g, b, a;
};
typedef std::vector<Pixel> Image;

inline float step(float edge, float x) {
    return x < edge ? 0.0f : 1.0f;
}
inline float min(float x, float y) {
    return std::min(x, y);
}
inline float max(float x, float y) {
    return std::max(x, y);
}
inline float clamp(float x, float minVal, float maxVal) {
    return min(max(x, minVal), maxVal);
}
struct vec2 {
    vec2() {}
    explicit vec2(float s) : x(s), y(s) {}
    vec2(float x, float y) : x(x), y(y) {}

    float operator[](int i) const { return *(&x + i); }
    float &operator[](int i)  { return *(&x + i); }
    vec2 operator-() { return vec2(-x, -y); }
    vec2 &operator+=(const vec2 &rhs) { x += rhs.x; y += rhs.y; return *this; }
    vec2 &operator-=(const vec2 &rhs) { x -= rhs.x; y -= rhs.y; return *this; }
    vec2 &operator*=(float rhs) { x *= rhs; y *= rhs; return *this; }

    float x, y;
};
inline vec2 operator+(float x, const vec2 &y) {
    return vec2(x + y[0], x + y[1]);
}
inline vec2 operator+(const vec2 &x, float y) {
    return vec2(x[0] + y, x[1] + y);
}
inline vec2 operator+(const vec2 &x, const vec2 &y) {
    return vec2(x[0] + y[0], x[1] + y[1]);
}
inline vec2 operator-(float x, const vec2 &y) {
    return vec2(x - y[0], x - y[1]);
}
inline vec2 operator-(const vec2 &x, float y) {
    return vec2(x[0] - y, x[1] - y);
}
inline vec2 operator-(const vec2 &x, const vec2 &y) {
    return vec2(x[0] - y[0], x[1] - y[1]);
}
inline vec2 operator*(float x, const vec2 &y) {
    return vec2(x * y[0], x * y[1]);
}
inline vec2 operator*(const vec2 &x, float y) {
    return vec2(x[0] * y, x[1] * y);
}
inline vec2 operator*(const vec2 &x, const vec2 &y) {
    return vec2(x[0] * y[0], x[1] * y[1]);
}
inline float dot(const vec2 &x, const vec2 &y) {
    return x[0]*y[0] + x[1]*y[1];
}
inline float length(const vec2 &x) {
    return std::sqrt(dot(x, x));
}
inline vec2 normalize(const vec2 &x) {
    float f = 1.0f / length(x);
    return f * x;
}
inline vec2 abs(const vec2 &x) {
    return vec2(x[0] < 0.0f ? -x[0] : x[0], x[1] < 0.0f ? -x[1] : x[1]);
}
inline vec2 min(const vec2 &x, float y) {
    return vec2(min(x[0], y), min(x[1], y));
}
inline vec2 min(const vec2 &x, const vec2 &y) {
    return vec2(min(x[0], y[0]), min(x[1], y[1]));
}
inline vec2 max(const vec2 &x, float y) {
    return vec2(max(x[0], y), max(x[1], y));
}
inline vec2 max(const vec2 &x, const vec2 &y) {
    return vec2(max(x[0], y[0]), max(x[1], y[1]));
}
inline vec2 floor(const vec2 &x) {
    return vec2(std::floor(x[0]), std::floor(x[1]));
}
inline vec2 step(float edge, const vec2 &x) {
    return vec2(step(edge, x[0]), step(edge, x[1]));
}
inline vec2 step(const vec2 &edge, const vec2 &x) {
    return vec2(step(edge[0], x[0]), step(edge[1], x[1]));
}
struct vec3 {
    vec3() {}
    explicit vec3(float s) : x(s), y(s), z(s) {}
    vec3(float x, float y, float z) : x(x), y(y), z(z) {}

    float operator[](int i) const { return *(&x + i); }
    float &operator[](int i)  { return *(&x + i); }
    vec3 operator-() { return vec3(-x, -y, -z); }
    vec3 &operator+=(const vec3 &rhs) { x += rhs.x; y += rhs.y; z += rhs.z; return *this; }
    vec3 &operator-=(const vec3 &rhs) { x -= rhs.x; y -= rhs.y; z -= rhs.z; return *this; }
    vec3 &operator*=(float rhs) { x *= rhs; y *= rhs; z *= rhs; return *this; }

    float x, y, z;
};
inline vec3 operator+(float x, const vec3 &y) {
    return vec3(x + y[0], x + y[1], x + y[2]);
}
inline vec3 operator+(const vec3 &x, float y) {
    return vec3(x[0] + y, x[1] + y, x[2] + y);
}
inline vec3 operator+(const vec3 &x, const vec3 &y) {
    return vec3(x[0] + y[0], x[1] + y[1], x[2] + y[2]);
}
inline vec3 operator-(float x, const vec3 &y) {
    return vec3(x - y[0], x - y[1], x - y[2]);
}
inline vec3 operator-(const vec3 &x, float y) {
    return vec3(x[0] - y, x[1] - y, x[2] - y);
}
inline vec3 operator-(const vec3 &x, const vec3 &y) {
    return vec3(x[0] - y[0], x[1] - y[1], x[2] - y[2]);
}
inline vec3 operator*(float x, const vec3 &y) {
    return vec3(x * y[0], x * y[1], x * y[2]);
}
inline vec3 operator*(const vec3 &x, float y) {
    return vec3(x[0] * y, x[1] * y, x[2] * y);
}
inline vec3 operator*(const vec3 &x, const vec3 &y) {
    return vec3(x[0] * y[0], x[1] * y[1], x[2] * y[2]);
}
inline float dot(const vec3 &x, const vec3 &y) {
    return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
}
inline float length(const vec3 &x) {
    return std::sqrt(dot(x, x));
}
inline vec3 normalize(const vec3 &x) {
    float f = 1.0f / length(x);
    return f * x;
}
inline vec3 cross(const vec3 &x, const vec3 &y) {
    return vec3(x[1]*y[2] - y[1]*x[2], x[2]*y[0] - y[2]*x[0], x[0]*y[1] - y[0]*x[1]);
}
inline vec3 abs(const vec3 &x) {
    return vec3(x[0] < 0.0f ? -x[0] : x[0], x[1] < 0.0f ? -x[1] : x[1], x[2] < 0.0f ? -x[2] : x[2]);
}
inline vec3 min(const vec3 &x, float y) {
    return vec3(min(x[0], y), min(x[1], y), min(x[2], y));
}
inline vec3 min(const vec3 &x, const vec3 &y) {
    return vec3(min(x[0], y[0]), min(x[1], y[1]), min(x[2], y[2]));
}
inline vec3 max(const vec3 &x, float y) {
    return vec3(max(x[0], y), max(x[1], y), max(x[2], y));
}
inline vec3 max(const vec3 &x, const vec3 &y) {
    return vec3(max(x[0], y[0]), max(x[1], y[1]), max(x[2], y[2]));
}
inline vec3 floor(const vec3 &x) {
    return vec3(std::floor(x[0]), std::floor(x[1]), std::floor(x[2]));
}
inline vec3 step(float edge, const vec3 &x) {
    return vec3(step(edge, x[0]), step(edge, x[1]), step(edge, x[2]));
}
inline vec3 step(const vec3 &edge, const vec3 &x) {
    return vec3(step(edge[0], x[0]), step(edge[1], x[1]), step(edge[2], x[2]));
}
struct vec4 {
    vec4() {}
    explicit vec4(float s) : x(s), y(s), z(s), w(s) {}
    vec4(float x, float y, float z, float w) : x(x), y(y), z(z), w(w) {}

    float operator[](int i) const { return *(&x + i); }
    float &operator[](int i)  { return *(&x + i); }
    vec4 operator-() { return vec4(-x, -y, -z, -w); }
    vec4 &operator+=(const vec4 &rhs) { x += rhs.x; y += rhs.y; z += rhs.z; w += rhs.w; return *this; }
    vec4 &operator-=(const vec4 &rhs) { x -= rhs.x; y -= rhs.y; z -= rhs.z; w -= rhs.w; return *this; }
    vec4 &operator*=(float rhs) { x *= rhs; y *= rhs; z *= rhs; w *= rhs; return *this; }

    float x, y, z, w;
};
inline vec4 operator+(float x, const vec4 &y) {
    return vec4(x + y[0], x + y[1], x + y[2], x + y[3]);
}
inline vec4 operator+(const vec4 &x, float y) {
    return vec4(x[0] + y, x[1] + y, x[2] + y, x[3] + y);
}
inline vec4 operator+(const vec4 &x, const vec4 &y) {
    return vec4(x[0] + y[0], x[1] + y[1], x[2] + y[2], x[3] + y[3]);
}
inline vec4 operator-(float x, const vec4 &y) {
    return vec4(x - y[0], x - y[1], x - y[2], x - y[3]);
}
inline vec4 operator-(const vec4 &x, float y) {
    return vec4(x[0] - y, x[1] - y, x[2] - y, x[3] - y);
}
inline vec4 operator-(const vec4 &x, const vec4 &y) {
    return vec4(x[0] - y[0], x[1] - y[1], x[2] - y[2], x[3] - y[3]);
}
inline vec4 operator*(float x, const vec4 &y) {
    return vec4(x * y[0], x * y[1], x * y[2], x * y[3]);
}
inline vec4 operator*(const vec4 &x, float y) {
    return vec4(x[0] * y, x[1] * y, x[2] * y, x[3] * y);
}
inline vec4 operator*(const vec4 &x, const vec4 &y) {
    return vec4(x[0] * y[0], x[1] * y[1], x[2] * y[2], x[3] * y[3]);
}
inline float dot(const vec4 &x, const vec4 &y) {
    return x[0]*y[0] + x[1]*y[1] + x[2]*y[2] + x[3]*y[3];
}
inline float length(const vec4 &x) {
    return std::sqrt(dot(x, x));
}
inline vec4 normalize(const vec4 &x) {
    float f = 1.0f / length(x);
    return f * x;
}
inline vec4 abs(const vec4 &x) {
    return vec4(x[0] < 0.0f ? -x[0] : x[0], x[1] < 0.0f ? -x[1] : x[1], x[2] < 0.0f ? -x[2] : x[2], x[3] < 0.0f ? -x[3] : x[3]);
}
inline vec4 min(const vec4 &x, float y) {
    return vec4(min(x[0], y), min(x[1], y), min(x[2], y), min(x[3], y));
}
inline vec4 min(const vec4 &x, const vec4 &y) {
    return vec4(min(x[0], y[0]), min(x[1], y[1]), min(x[2], y[2]), min(x[3], y[3]));
}
inline vec4 max(const vec4 &x, float y) {
    return vec4(max(x[0], y), max(x[1], y), max(x[2], y), max(x[3], y));
}
inline vec4 max(const vec4 &x, const vec4 &y) {
    return vec4(max(x[0], y[0]), max(x[1], y[1]), max(x[2], y[2]), max(x[3], y[3]));
}
inline vec4 floor(const vec4 &x) {
    return vec4(std::floor(x[0]), std::floor(x[1]), std::floor(x[2]), std::floor(x[3]));
}
inline vec4 step(float edge, const vec4 &x) {
    return vec4(step(edge, x[0]), step(edge, x[1]), step(edge, x[2]), step(edge, x[3]));
}
inline vec4 step(const vec4 &edge, const vec4 &x) {
    return vec4(step(edge[0], x[0]), step(edge[1], x[1]), step(edge[2], x[2]), step(edge[3], x[3]));
}
//
// Description : Array and textureless GLSL 2D/3D/4D simplex
// noise functions.
// Author : Ian McEwan, Ashima Arts.
// Maintainer : ijm
// Lastmod : 20110822 (ijm)
// License : Copyright (C) 2011 Ashima Arts. All rights reserved.
// Distributed under the MIT License. See LICENSE file.
// https://github.com/ashima/webgl-noise
//
inline vec3 mod289(const vec3 &x) {
    return x - floor(x * (1.0f / 289.0f)) * 289.0f;
}
inline vec4 mod289(const vec4 &x) {
    return x - floor(x * (1.0f / 289.0f)) * 289.0f;
}
inline vec4 permute(const vec4 &x) {
    return mod289(((x * 34.0f) + 1.0f) * x);
}
inline vec4 taylorInvSqrt(const vec4 &r) {
    return 1.79284291400159f - 0.85373472095314f * r;
}
float snoise(const vec3 &v) {
    const vec2 C = vec2(1.0f/6.0f, 1.0f/3.0f);
    const vec4 D = vec4(0.0f, 0.5f, 1.0f, 2.0f);
    vec3 i = floor(v + dot(v, vec3(C.y)));
    vec3 x0 = v - i + dot(i, vec3(C.x));
    vec3 g = step(vec3(x0.y, x0.z, x0.x), x0);
    vec3 l = 1.0f - g;
    vec3 i1 = min(g, vec3(l.z, l.x, l.y));
    vec3 i2 = max(g, vec3(l.z, l.x, l.y));
    vec3 x1 = x0 - i1 + vec3(C.x);
    vec3 x2 = x0 - i2 + vec3(C.y);
    vec3 x3 = x0 - vec3(D.y);
    i = mod289(i);
    vec4 p = permute(permute(permute(i.z + vec4(0.0f, i1.z, i2.z, 1.0f)) + i.y + vec4(0.0f, i1.y, i2.y, 1.0f)) + i.x + vec4(0.0f, i1.x, i2.x, 1.0f));
    float n_ = 0.142857142857f;
    vec3 ns = n_ * vec3(D.w, D.y, D.z) - vec3(D.x, D.z, D.x);
    vec4 j = p - 49.0f * floor(p * ns.z * ns.z);
    vec4 x_ = floor(j * ns.z);
    vec4 y_ = floor(j - 7.0f * x_);
    vec4 x = x_ * ns.x + vec4(ns.y);
    vec4 y = y_ * ns.x + vec4(ns.y);
    vec4 h = 1.0f - abs(x) - abs(y);
    vec4 b0 = vec4(x.x, x.y, y.x, y.y);
    vec4 b1 = vec4(x.z, x.w, y.z, y.w);
    vec4 s0 = floor(b0)*2.0f + 1.0f;
    vec4 s1 = floor(b1)*2.0f + 1.0f;
    vec4 sh = -step(h, vec4(0.0f));
    vec4 a0 = vec4(b0.x, b0.z, b0.y, b0.w) + vec4(s0.x, s0.z, s0.y, s0.w)*vec4(sh.x, sh.x, sh.y, sh.y);
    vec4 a1 = vec4(b1.x, b1.z, b1.y, b1.w) + vec4(s1.x, s1.z, s1.y, s1.w)*vec4(sh.z, sh.z, sh.w, sh.w);
    vec3 p0 = vec3(a0.x, a0.y, h.x);
    vec3 p1 = vec3(a0.z, a0.w, h.y);
    vec3 p2 = vec3(a1.x, a1.y, h.z);
    vec3 p3 = vec3(a1.z, a1.w, h.w);
    vec4 norm = taylorInvSqrt(vec4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)));
    p0 *= norm.x;
    p1 *= norm.y;
    p2 *= norm.z;
    p3 *= norm.w;
    vec4 m = max(0.6f - vec4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0f);
    m = m*m;
    return 42.0f * dot(m*m, vec4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)));
}
inline float fbm(const vec3 &p) {
    return
        0.5000f * snoise(1.0f * p) +
        0.2500f * snoise(2.0f * p) +
        0.1250f * snoise(4.0f * p) +
        0.0625f * snoise(8.0f * p);
}
inline float box(const vec3 &p, const vec3 &b, float r) {
    return length(max(abs(p) - b, 0.0f)) - r;
}
inline float sphere(const vec3 &p, const vec3 &c, float r) {
    return length(p - c) - r;
}
float map(const vec3 &p, int &matID) {
    float dist, minDist = 100.0f;
    dist = box(p, vec3(1.5f, 1.0f, 1.5f), 0.03f);
    if (dist < minDist) {
        minDist = dist;
        matID = 1;
    }
    dist = sphere(p, vec3(0.0f, 1.0f, 0.0f), 1.0f);
    if (dist < minDist) {
        minDist = dist;
        matID = 2;
    }
    dist = dot(p, vec3(0.0f, 1.0f, 0.0f));
    if (dist < minDist) {
        minDist = dist;
        matID = 3;
    }
    return minDist;
}
float is(const vec3 &ro, const vec3 &rd, int &matID) {
    for (float t = 0.0f; t < 40.0f;) {
        float d = map(ro + rd*t, matID);
        if (d < 0.001f)
            return t;
        t += d;
    }
    matID = -1;
    return -1.0f;
}
float shadow(const vec3 &ro, const vec3 &rd) {
    float r = 1.0f;
    int m;
    for (float t = 0.01f; t < 10.0f;) {
        float d = map(ro + rd*t, m);
        if (d < 0.0005f)
            return 0.0f;
        r = min(r, 16.0f*d/t);
        t += d;
    }
    return r;
}
vec3 computeNormal(const vec3 &p) {
    int m;
    vec3 n;
    n.x = map(vec3(p.x + 0.001f, p.y, p.z), m) - map(vec3(p.x - 0.001f, p.y, p.z), m);
    n.y = map(vec3(p.x, p.y + 0.001f, p.z), m) - map(vec3(p.x, p.y - 0.001f, p.z), m);
    n.z = map(vec3(p.x, p.y, p.z + 0.001f), m) - map(vec3(p.x, p.y, p.z - 0.001f), m);
    return n;
}
vec3 computeBump(const vec3 &p) {
    vec3 b;
    float d = fbm(p);
    b.x = fbm(vec3(p.x + 0.001f, p.y, p.z)) - d;
    b.y = fbm(vec3(p.x, p.y + 0.001f, p.z)) - d;
    b.z = fbm(vec3(p.x, p.y, p.z + 0.001f)) - d;
    return b;
}

void computeColor(Pixel &out, float x, float y) {
    vec3 ro = vec3(2.0f, 3.0f, 3.0f);
    vec3 az = normalize(vec3(0.0f) - ro);
    vec3 ax = normalize(cross(vec3(0.0f, 1.0f, 0.0f), az));
    vec3 ay = normalize(cross(az, ax));
    vec3 rd = normalize(ax*x + ay*y + az*1.5f);

    vec3 c = vec3(0.0f);
    int matID;
    float d = is(ro, rd, matID);
    if (d > 0.0f) {
        vec3 p = ro + d*rd;
        vec3 n = computeNormal(p);
        /*
        if (matID == 2) {
            vec3 b = computeBump(8.0f * p);
            n = normalize(n + 0.5f * b);
        }
        else {*/
            n = normalize(n);
        //}
        //c = 0.5f*(1.0f+n);

            /*
        // AO
        float ao;
        float totao = 0.0f;
        float sca = 10.0f;
        for (int aoi = 0; aoi < 5; ++aoi) {
            float hr = 0.01f + 0.015f * static_cast<float>(aoi * aoi);
            vec3 aopos = n * hr + p;
            int m;
            float dd = map(aopos, m);
            ao = -(dd - hr);
            totao += ao*sca;
            sca *= 0.5f;
        }
        ao = 1.0f - clamp(totao, 0.0f, 1.0f);
        */

        vec3 rgb = vec3(1.0f);

        /*
        if (matID == 2) {
            float fb = fbm(8.0f * p);
            rgb = vec3(0.64f, 0.61f, 0.59f) + fb*vec3(0.21f, 0.19f, 0.19f);
        }
        */

        vec3 lightDir0 = normalize(vec3(10.0f, 8.0f, -6.0f) - p);
        vec3 lightDir1 = normalize(vec3(-12.0f, 19.0f, 6.0f) - p);

        c += 0.7f*max(0.4f + 0.6f*dot(n, lightDir0)*shadow(p, lightDir0), 0.0f)*rgb;
        c += 0.3f*max(0.4f + 0.6f*dot(n, lightDir1)*shadow(p, lightDir1), 0.0f)*rgb;
        //c *= ao;
    }

    out.r = c.x;
    out.g = c.y;
    out.b = c.z;
    out.a = 1.0f;
}

int main(int argc, char *argv[]) {
    Image image(IMAGE_WIDTH * IMAGE_HEIGHT, Pixel(1.0f, 0.0f, 0.0f, 1.0f));

    for (int j = 0; j < IMAGE_HEIGHT; ++j) {
        for (int i = 0; i < IMAGE_WIDTH; ++i) {
            float x = (i + 0.5f) / IMAGE_WIDTH;
            x = (x - 0.5f) * 2.0f * 1.777f;
            float y = (j + 0.5f) / IMAGE_HEIGHT;
            y = (y - 0.5f) * 2.0f;
            computeColor(image[j * IMAGE_WIDTH + i], x, y);
        }
    }

    std::ofstream file("test0.tga", std::ios::binary);

    const char tgaHeader[18] = { 0,0,2,0,0,0,0,0,0,0,0,0,IMAGE_WIDTH & 0x00ff,(IMAGE_WIDTH & 0xff00) >> 8,IMAGE_HEIGHT & 0x00ff,(IMAGE_HEIGHT & 0xff00) >> 8,32,0 };
    file.write(tgaHeader, sizeof(tgaHeader));

    for (Image::iterator i = image.begin(); i != image.end(); ++i) {
        char rgba[4];
        rgba[0] = static_cast<char>(255.0f * i->b);
        rgba[1] = static_cast<char>(255.0f * i->g);
        rgba[2] = static_cast<char>(255.0f * i->r);
        rgba[3] = static_cast<char>(255.0f * i->a);
        file.write(rgba, sizeof(rgba));
    }

    return 0;
}
